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ABSTRACT 

This  is  a  review  of  molecular  modelling  techniques  which  may  be  applied  to  studies  of 
energetic  materials.  It  focusses  on  ab  initio  ('first-principles')  molecular  orbital 
calculations,  since  these  methods  offer  the  greatest  accuracy.  Since  ab  initio  calculations 
are  very  computer-intensive,  approximate  MO  methods  are  also  discussed,  which  offer 
reasonably  accurate  predictions  with  reduced  calculation  times.  These  approximate 
methods  include  density  functional  theory  and  'layered'  techniques  which  combine 
different  levels  of  theoretical  sophistication  into  one  calculation. 
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Overview  of  Molecular  Modelling  and 
Ab  initio  Molecular  Orbital  Methods 
Suitable  for  Use  with  Energetic  Materials 


Executive  Summary 

Qualifying  energetic  materials  (EM)  and  formulations  for  military  use  is  a  tedious, 
expensive  and  hazardous  process.  It  may  be  possible  to  replace  much  of  the  redundant 
empirical  testing  of  EM,  either  in  service  or  under  consideration,  by  computer 
simulation.  Therefore,  accurate  theoretical  models  of  EM  performance  and  stability  can 

•  Accelerate  research  and  development  of  new  materials  with  unproved 

performance  and  handling  characteristics 

•  Reduce  maintenance  costs  for  energetic  materials  in  service  and 

•  Reduce  risk  to  personnel  testing,  handling  or  using  energetic  materials. 

Many  relevant  characteristics  of  energetic  materials  -  for  instance,  the  performance  and 
sensitivity  of  explosives,  or  the  stability  and  signatures  of  propellants  -  are  controlled 
by  their  chemical  properties.  These  characteristics  depend  upon  properties  such  as 
molecular  structures  and  reaction  energies  that  can  be  predicted  with  computational 
chemistry  techniques  which  focus  on  molecular  modelling. 

All  molecular  modelling  techniques  can  be  classified  under  three  general  categories:  ab 
initio  electronic  structure  calculations,  semi-empirical  methods,  and  molecular 
mechanics.  Ab  initio  or  'first  principles'  electronic  structure  methods  are  based  upon 
quantum  mechanics  and  therefore  provide  the  most  accurate  and  consistent 
predictions  for  chemical  systems.  However  ab  initio  methods  are  extremely  computer¬ 
intensive.  Semi-empirical  methods  are  also  founded  upon  quantum  mechanics,  but 
speed  computation  by  replacing  some  explicit  calculations  with  approximations  based 
upon  experimental  data.  Molecular  mechanics  techniques  are  purely  empirical 
methods  based  upon  the  principles  of  classical  physics,  and  as  such  are 
computationally  fast.  Molecular  mechanics  methods  completely  neglect  explicit 
treatment  of  electronic  structure,  and  are  therefore  severely  limited  in  scope;  however, 
they  often  provide  the  only  means  with  which  to  study  very  large  chemical  systems 
(e.g.,  polymers  or  solutions)  or  non-homogeneous  mixtures  like  those  typically  used  in 
conventional  explosives,  pyrotechnics  or  propellants  formulations. 

Ab  initio  methods  are  capable  of  high  accuracy  predictions  over  a  wide  range  of 
systems.  Rapid  advances  in  computer  technology  are  making  ab  initio  methods 
increasingly  more  practical  for  use  with  realistic  chemical  systems.  Similarly, 
computationally  'cheaper'  techniques  such  as  density  functional  calculations  and 
'layered'  methods  such  as  ONIOM,  are  continually  being  refined.  These  methods  show 


promise  of  providing  consistent  and  accurate  chemical  predictions  for  complicated 
systems  requiring  explicit  treatment  of  electronic  structure,  such  as  energetic  molecules 
that  contain  azides  or  nitrogen  oxides. 

This  document  describes  the  theoretical  basis  of  molecular  modelling  techniques  which 
may  be  applied  to  energetic  materials.  A  subsequent  publication  will  describe  in  detail 
the  application  of  electronic  structure  methods  to  the  study  of  energetic  materials. 
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1.  Introduction 


Computational  chemistry  may  be  defined  as  the  application  of  mathematical  and 
theoretical  principles  to  the  solution  of  chemical  problems.  Molecular  modelling,  a 
subset  of  computational  chemistry,  concentrates  on  predicting  the  behaviour  of 
individual  molecules  within  a  chemical  system.  The  most  accurate  molecular  models 
use  ab  initio  or  'first  principles'  electronic  structure  methods,  based  upon  the  principles 
of  quantum  mechanics,  and  are  generally  very  computer-intensive.  However,  due  to 
advances  in  computer  storage  capacity  and  processor  performance,  molecular 
modelling  has  been  a  rapidly  evolving  and  expanding  field,  to  the  point  that  it  is  now 
possible  to  solve  relevant  problems  in  an  acceptable  amount  of  time. 

The  types  of  predictions  possible  for  molecules  and  reactions  include  [1]: 

•  Heats  of  formation 

•  Bond  and  reaction  energies 

•  Molecular  energies  and  structures  (thermochemical  stability) 

•  Energies  and  structures  of  transition  states  (activation  energies) 

•  Reaction  pathways,  kinetics  and  mechanisms 

•  Charge  distribution  in  molecules  (reactive  sites) 

•  Substituent  effects 

•  Electron  affinities  and  ionisation  potentials 

•  Vibrational  frequencies  (IR  and  Raman  spectra) 

•  Electronic  transitions  (UV/ Visible  spectra) 

•  Magnetic  shielding  effects  (NMR  spectra) 

Prediction  of  these  properties  has  many  applications  in  energetic  materials  research, 
including  studies  of  synthesis  pathways,  reaction  products,  initiation  mechanisms,  and 
exhaust  plume  signature.  This  report  presents  a  review  of  computational  chemistry 
techniques,  focussing  on  electronic  structure  methods.  Electronic  structure  methods, 
particularly  ab  initio  calculations,  are  capable  of  consistent  predictions  with  high 
accuracy  (±  20  kj/mol)1  over  a  wide  range  of  systems  -  a  critical  prerequisite  for  the 
successful  modelling  of  energetic  materials. 

An  example  of  the  application  of  molecular  modelling  to  energetic  materials  is 
calculation  of  the  detonation  properties  of  a  new  explosive  FOX-7  (l,l-diamino-2,2- 
dinitroethylene).  The  ideal  detonation  properties  of  any  explosive  formulation  can  be 
calculated  from  the  principles  of  chemical  thermodynamics  using  only  the  heats  of 
formation  of  the  ingredients  as  input.  The  experimental  heat  of  formation  of  FOX-7  is 


1  The  unit  'kj/mol'  is  a  convenient  scale  for  discussing  molecular  energies.  For  instance,  bonds 
within  molecules  typically  have  energies  of  a  few  hundred  kj/mol;  bonds  between  molecules 
(e.g.,  hydrogen  bonds)  have  energies  of  tens  of  kj/mol;  and  bonds  between  atoms  of  inert  gases 
(e.g.,  He,  Ne,  Ar)  have  energies  no  more  than  a  few  kj/  mol. 
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-134  kj/mol  [2],  and  the  calculated  heat  of  formation  using  the  PM-3  semi-empirical 
method  is  -53  kj/mol.  Estimates  of  the  velocity  of  detonation  (VoD)  and  Chapman- 
Jouguet  detonation  pressure  (Pc-j)  for  FOX-7  calculated  with  the  CHEETAH 
thermochemical  code  [3]  are  presented  in  Table  1.1.  In  this  case,  it  is  clear  that  the  PM-3 
heat  of  formation  is  accurate  enough  to  yield  values  of  VoD  and  Pc-j  that  are  very  close 
to  those  estimated  from  the  experimental  heat  of  formation.  Heats  of  formation 
calculated  with  ab  initio  methods  are  nearly  identical  to  the  experimental  value. 

Table  1.1  FOX-7  detonation  parameters  estimated  with  CHEETAH. 


AHf 

CHEETAH  v(1.39) 

CHEETAH  vl  .41 

VoD 

Pc-j 

VoD 

Pc-J 

-53  kj/mol 

9130  m/s 

37.1  GPa 

9040  m/s 

36.0  GPa 

-134  kj/mol 

9090  m/s 

36.6  GPa 

8870  m/  s 

34.0  GPa 

Ab  initio ,  semi-empirical  and  molecular  mechanics  methods  have  been  used  extensively 
to  study  the  chemistry  of  energetic  materials  in  explosives,  propellants  and 
pyrotechnics.  A  comprehensive  review  of  published  studies  is  beyond  the  scope  of  this 
paper,  but  will  be  summarised  in  a  subsequent  technical  report  [4]. 


2.  Survey  of  Computational  Chemistry  Methods 

All  molecular  modelling  techniques  can  be  classified  under  three  major  categories: 
ab  initio  electronic  structure  calculations,  semi-empirical  methods  and  molecular 
mechanics.  General  characteristics  for  each  method  are  summarised  in  Table  2.1. 


2.1  Ab  initio  electronic  structure  methods 

Of  the  three,  ab  initio  molecular  orbital  methods  are  the  most  accurate  and  consistent 
because  they  provide  the  best  mathematical  approximation  to  the  actual  system.  The 
term  ab  initio  implies  that  the  computations  are  based  solely  on  the  laws  of  quantum 
mechanics,  the  masses  and  charges  of  electrons  and  atomic  nuclei,  and  the  values  of 
fundamental  physical  constants,  such  as  the  speed  of  light  (c  =  2.998x10s  m/s)  or 
Planck's  constant  (h  =  6.626X1034  J-s),  and  contain  no  approximations.  'Molecular 
orbital'  methods  solve  Schrodinger's  equation  for  the  chemical  system  using  a  'basis 
set'  of  functions  that  satisfy  a  series  of  rigorous  mathematical  approximations.  Ab  initio 
concepts  are  outlined  in  Chapter  3. 

Ab  initio  molecular  orbital  calculations  are  specified  by  a  'model  chemistry.'  The  model 
chemistry  includes  the  choice  of  method  and  basis  set,  the  general  structure  and 
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electronic  state  of  the  molecular  system  under  study  (e.g.,  charge  and  spin  states),  and 
the  treatment  of  electron  spin2.  Molecular  properties  can  be  assessed  from  a  user- 
specified  input  (single-point  energy  calculation  or  SPE),  or  the  molecule  can  be  allowed 
to  relax  to  a  minimum  energy  configuration  (geometry  optimisation). 

Table  2.1.  Synopsis  of  molecular  modelling  techniques. 


Method 

Advantages 

Disadvantages 

Best  for 

Ab  initio 

•  Useful  for  a 

•  Computationally 

•  Small  systems 

•  Uses  quantum 

broad  range  of 

expensive 

(tens  of  atoms) 

physics 

systems 

•  Electronic 

•  Mathematically 

•  Does  not  depend 

transitions 

rigorous:  no 

on  experimental 

•  Systems  without 

empirical 

data 

experimental 

parameters 

•  Calculates 

data 

transition  states 

•  Systems 

and  excited  states 

requiring  high 

accuracy. 

Semi-empirical 

•  Uses  quantum 

•  Less  demanding 

•  Requires  ab  initio 

•  Medium-sized 

physics 

computationally 

or  experimental 

systems 

•  Uses 

than  ab  initio 

data  for 

(hundreds  of 

experimental 

methods. 

parameters. 

atoms). 

parameters 

•  Calculates 

•  Less  rigorous  than 

•  Electronic 

•  Uses  extensive 

transition  states 

ab  initio  methods. 

transitions. 

approximations 

and  excited 

states. 

Molecular 

Mechanics 

•  Uses  classical 

•  Computationally 

•  Does  not  calculate 

•  Large  systems 

physics 

'cheap:'  fast  and 

electronic 

(thousands  of 

•  Relies  on  force 

useful  with 

properties 

atoms) 

field  with 

limited  computer 

•  Requires  ab  initio 

•  Systems  or 

embedded 

resources. 

or  experimental 

processes  that  do 

empirical 

•  Can  be  used  for 

data  for 

not  involve  bond 

parameters. 

large  molecules 

parameters 

breaking. 

like  enzymes. 

•  Commercial 

software 

applicable  to  a 

limited  range  of 

molecules 

Ab  initio  molecular  orbital  computations  can  provide  accurate  quantitative  predictions 
of  chemical  properties  for  a  wide  range  of  molecular  systems.  However,  they  place  a 
considerable  demand  on  computer  resources.  The  choice  of  theoretical  method  and 
basis  set  determine  the  duration  of  the  calculation;  thus,  a  sophisticated  method  and  a 


2  Electrons  may  be  assigned  singly  (unrestricted)  or  in  pairs  of  opposite  spin  (restricted)  to  the 
molecular  orbitals  that  make  up  the  total  wavefunction  of  the  system. 
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large  basis  set  will  provide  more  accurate  results,  but  will  also  require  more  computer 
resources.  Hence  timely  calculations  on  large  molecules  (n  ~  10-30  atoms3)  may  only  be 
possible  using  the  basic  theory  ( i.e .,  Hartree  Fock  approximation)  with  minimal  basis 
sets,  whereas  calculations  on  chemical  reactions  between  simple  diatomic  molecules 
can  be  performed  with  state-of  the-art  model  chemistries.  For  very  large  systems  («  ~ 
50  atoms  or  more),  timely  results  are  only  possible  with  less  sophisticated  models,  like 
those  described  in  sections  2.2  and  2.3  below. 


2.2  Semi-empirical  methods 

Semi-empirical  methods  increase  the  speed  of  computation  by  using  approximations  of 
ab  initio  techniques  (e.g.,  by  limiting  choices  of  molecular  orbitals  or  considering  only 
valence  electrons)  which  have  been  fitted  to  experimental  data  (for  instance,  structures 
and  formation  energies  of  organic  molecules).  Until  recently,  the  size  of  many  energetic 
molecules  placed  them  beyond  the  scope  of  ab  initio  calculations,  so  preliminary 
theoretical  studies  were  performed  using  semi-empirical  techniques  [3-11] .  However, 
semi-empirical  methods  have  been  calibrated  to  typical  organic  or  biological  systems 
and  tend  to  be  inaccurate  for  problems  involving  hydrogen-bonding,  chemical 
transitions  or  nitrated  compounds  [5,7,12,13]. 

Several  semi-empirical  methods  are  available  and  appear  in  commercially  available 
computational  chemistry  software  packages  such  as  HyperChem  [14]  and  Chem3D 
[15].  Some  of  the  more  common  semi-empirical  methods  can  be  grouped  according  to 
their  treatment  of  electron-electron  interactions  [14]. 

2.2.1  The  extended  Hiickel  method 

Extended  Hiickel  calculations  neglect  all  electron-electron  interactions,  making  them 
computationally  fast  but  not  very  accurate.  The  model  provides  a  qualitative  estimate 
of  the  shapes  and  relative  energies  of  molecular  orbitals,  and  approximates  the  spatial 
distribution  of  electron  density.  Extended  Hiickel  models  are  good  for  chemical 
visualisation  and  can  be  applied  to  'frontier  orbital'  treatments  of  chemical  reactivity. 

2.2.2  Neglect  of  differential  overlap  (NDO) 

NDO  models  neglect  some  but  not  all  of  the  electron-electron  interactions.  The  Hartree- 
Fock  Self-Consistent  Field  (HF-SCF)  method  (see  section  3.2)  is  used  to  solve  the 
Schrodinger  equation  with  various  approximations: 

•  Complete  NDO  (CNDO)  -  the  product  of  two  atomic  orbitals  on  different  atoms  is 
set  equal  to  zero  everywhere.  In  this  case,  repulsion  between  electrons  in 
different  orbitals  depends  only  upon  the  nature  of  the  atoms  involved,  and  not 


3  Obviously,  n  is  constantly  increasing  as  computing  power  is  improved. 
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on  the  particular  orbitals.  Because  CNDO  neglects  almost  all  description  of 
electron  exchange  properties,  it  does  not  differentiate  between  states  that  have 
the  same  electronic  configuration,  but  different  values  of  electron  spin. 

•  Intermediate  NDO  (INDO)  -  differential  overlap  between  orbitals  on  the  same 
atom  are  taken  into  account  in  the  description  of  electron-electron  repulsion,  but 
differential  overlap  between  orbitals  on  different  atoms  is  neglected. 

•  Modified  INDO,  version  3  (MINDO/3)  -  reparameterised  version  of  INDO 
optimised  to  predict  good  enthalpies  of  formation  and  reasonable  molecular 
geometries  for  a  range  of  chemical  systems,  in  particular,  sulphur-containing 
compounds,  carbocations,  and  polynitro  organic  compounds  [16]. 

•  Zemer’s  INDO  methods  (ZINDO/1  and  ZINDO/S)  -  Michael  Zemer's  (University 
of  Florida)  versions  of  INDO  developed  for  use  with  molecular  systems 
containing  transition  metals.  ZINDO/1  is  optimised  to  predict  molecular 
geometries  and  ZINDO/S  is  optimised  to  predict  UV  spectra. 

2.2.3  Neglect  of  diatomic  differential  overlap  (NDDO) 

NDDO  methods  build  upon  the  INDO  model  by  including  the  overlap  density 
between  two  orbitals  on  one  atom  interacting  with  the  overlap  density  between  two 
orbitals  on  the  same  or  another  atom. 

•  Modified  NDO  (MNDO)  -  a  method  introduced  to  correct  some  of  the  problems 
associated  with  MLNDO/3  [14].  MNDO  does  not  work  well  for  sterically 
crowded  molecules,  four-membered  rings,  hydrogen  bonding,  hypervalent 
compounds,  nitro  groups  and  peroxides.  In  general,  MNDO  overestimates 
activation  barriers  to  chemical  reactions  [15]. 

•  Austin  Method,  version  1  (AMI)  -  a  reparameterised  version  of  MNDO  which 
includes  changes  in  nuclear  repulsion  terms  [15].  Although  more  accurate  than 
MNDO,  AMI  does  not  handle  phosphorus-oxygen  bonds,  nitro  compounds  and 
peroxide  bonds  [14]. 

•  Parameterisation  Model,  version  3  (PM3)  -  a  second  reparameterisation  of  MNDO, 
functionally  .similar  to  AMI,  but  with  some  significant  improvements.  PM3  is  a 
recently  developed  semi-empirical  method  that  may  contain  as  yet  undiscovered 
defects  [15]. 


2.3  Molecular  mechanics 

Molecular  mechanics  (MM)  is  often  die  only  feasible  means  with  which  to  model  very 
large  and  non-symmetric  chemical  systems  such  as  proteins  or  polymers.  Molecular 
mechanics  is  a  purely  empirical  method  that  neglects  explicit  treatment  of  electrons, 
relying  instead  upon  the  laws  of  classical  physics  to  predict  the  chemical  properties  of 
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moleaxles.  As  a  result,  MM  calculations  cannot  deal  with  problems  such  as  bond 
breaking  or  formation,  where  electronic  or  quantum  effects  dominate.  Furthermore, 
MM  models  are  wholly  system-dependent;  MM  energy  predictions  tend  to  be 
meaningless  as  absolute  quantities,  and  are  generally  useful  only  for  comparative 
studies.  Despite  these  shortcomings,  MM  bridges  the  gap  between  quantum  and 
continuum  mechanics,  and  has  been  used  quite  extensively  to  study  'mesoscopic' 
effects  in  energetic  materials.  Applications  include  modelling  reaction  and  dissociation 
on  classical  potential  energy  surfaces  [17-20],  studies  of  equilibrium  crystal  properties 
(e.g.,  density,  packing,  specific  heats)  [21-27],  dynamic  investigations  of  shock 
interactions  with  crystals  and  defects  [28,29],  and  simulating  detonation  in  molecular 
crystals  [30-33]. 

The  basic  assumptions  of  typical  molecular  mechanics  methods  are  listed  below. 

•  Each  atom  (i.e.,  electrons  and  nucleus)  is  represented  as  one  particle  with  a 
characteristic  mass. 

•  A  chemical  bond  is  represented  as  a  'spring,'  with  a  characteristic  force  constant 
determined  by  the  potential  energy  of  interaction  between  the  two  participating 
atoms.  Potential  energy  functions  can  describe  intramolecular  bond  stretching, 
bending  and  torsion,  or  infcrmolecular  phenomena  such  as  electrostatic 
interactions  or  van  der  Waals  forces. 

•  The  potential  energy  functions  rely  on  empirically  derived  parameters  obtained 
from  experiments  or  from  other  calculations. 


Current  molecular  mechanics  models  are  characterised  by  the  set  of  potential  energy 
functions  used  to  describe  the  chemical  forces.  These  force  fields  depend  upon: 

•  Atomic  displacements  (i.e.,  bond  lengths); 

•  Atom  types,  that  is,  the  characteristics  of  an  element  within  a  specific  chemical 
context  (e.g.,  a  carbonyl  carbon  versus  a  methyl  carbon);  and 

•  One  or  more  parameter  sets  relating  atom  types  and  bond  characteristics  to 
empirical  data. 


Examples  of  MM  force  fields  in  common  use  are: 

•  AMBER  (Assisted  Model  Building  with  Energy  Refinement)  -  primarily  designed  for 
the  study  of  biomolecules  such  as  proteins  and  nucleotides  [34]. 

•  CHARMM  (Chemistry  at  HARvard  Molecular  Mechanics)  -  primarily  designed  for 
biological  and  pharmaceutical  study,  but  has  also  been  applied  to  micelles  and 
self-assembling  macromolecules  [35]. 
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•  MMx  (MM2,  MM3,  etc)  -  optimised  for  structural  and  thermodynamic  studies  of 
small  non-polar  molecules  [14].  MMx  force  fields  include  third-  and  fourth-order 
corrections  to  standard  quadratic  fits  for  the  potential  energy  surfaces  of  bonds 
and  bond  angles,  thus  allowing  for  non-harmonic  effects  in  molecular 
vibrations.  The  various  MMx  versions  differ  primarily  in  their 
parameterisations.  The  higher  versions  tend  to  be  more  modem  and  address 
deficiencies  in  their  predecessors.  However,  for  the  newer  versions  such  as 
MM4,  parameters  may  not  be  available  for  all  classes  of  molecules. 

•  OPLS  (Optimised  Potentials  for  Liquid  Simulations)  -  optimised  for  reproducing  the 
physical  properties  of  biomolecules  in  liquid  solutions  [14]. 


The  authors  are  not  aware  of  any  published  study  discussing  the  use  of  AMBER, 
CHARMM,  MMx  or  OPLS  with  energetic  materials.  Since  these  packages  are  optimised 
for  biochemistry  and  pharmaceutical  applications,  it  is  unlikely  that  they  will 
accurately  reproduce  the  behaviour  of  energetic  materials  without  further 
modification.  However,  it  is  likely  they  can  be  use  for  limited  applications  with  only 
slight  modification. 

Basic  techniques  of  molecular  dynamics  are  described  in  texts  such  as  Rapaport  [36] 
and  Haile  [37],  and  a  general  review  of  molecular  mechanics  applications  in  several 
areas  of  chemistry  is  given  by  Rappe  [38]. 


2.4  Hybrid  methods 

It  is  worth  mentioning  here  recently  developed  hybrid  methods  which  combine  high- 
level  quantum  mechanical  calculations  on  a  small  part  of  a  system  with  a  lower-level 
method  on  the  rest  of  the  system.  Thus  for  large  clusters  or  macromolecules,  accurate 
calculations  can  be  carried  out  on  the  area  of  interest  without  ignoring  or  making 
unnecessary  assumptions  about  the  remainder  of  the  system 

Three  methods  have  been  developed  and  extensively  used  by  Morokuma  et  al.  at 
Emory  University  (Atlanta,  GA).  These  methods  and  their  applications  are  likely  to 
receive  increasing  attention  in  the  next  few  years: 

•  Integrated  Molecular  Orbital  +  Molecular  Mechanics  (IMOMM)  -  a  two-layer  method 
in  which  a  high  level  MO  calculation  is  combined  with  molecular  mechanics  [39]. 
For  example,  the  IMOMM  method  has  been  used  to  investigate  organometallic- 
catalysed  polymerisation  reactions  [40,41]. 

•  Integrated  Molecular  Orbital  +  Molecular  Orbital  (IMOMO)  -  a  two-layer  method 
which  combines  high  level  and  low-level  MO  calculations  [42,43].  This  method  can 
be  used  for  calculating  bond  dissociation  energies  of  large  molecules  [44]. 

•  "Our  own  N-layered  Integrated  molecular  Orbital  molecular  Mechanics  (ONIOM) 
method"  -  a  technique  that  has  largely  superseded  IMOMM  and  IMOMO,  which  are 
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in  fact,  subsets  of  ONIOM.  ONIOM  is  a  three-level  general  method,  allowing  the 
combination  of  any  high-level  method  with  any  low-  or  'medium' -level  method 
[45].  ONIOM  has  been  applied  to  calculations  of  bond  breaking  in  fullerenes  [46]. 


3.  Ab  initio  Molecular  Orbital  Theory 

What  follows  is  a  brief  description  of  the  fundamental  concepts  that  underlie  the 
algorithms  of  typical  off-the-shelf  ab  initio  molecular  orbital  software.  More 
comprehensive  discussion  of  molecular  orbital  theory  may  be  found  in  selected 
references  [12,13,16,47-49]. 


3.1  The  physico-chemical  model 

The  basis  of  electronic  structure  methods  is  the  assumption  that  all  chemistry  can  be 
described  in  terms  of  the  interactions  between  electronic  charges  within  molecules. 
Hence,  chemical  bonds  can  be  loosely  defined  as  a  redistribution  of  electronic  charge  that 
stabilises  the  molecule  with  respect  to  a  collection  of  its  (isolated)  constituent  atoms. 

Relative  stabilities  are  expressed  in  terms  of  the  total  energy  of  the  system,  which  is 
defined  by  a  differential  equation, 

H  =  T+ V,  (3.1.1) 

where  H  is  the  Hamiltonian  operator  representing  the  sum  of  kinetic  ( T )  and  potential 
( V )  energies.  In  quantum  mechanical  systems,  the  kinetic  energy  of  a  particle  is 

.  -h2  , 

T  =  — V2,  (3.1.2) 

2m 


where  m  is  the  mass  of  the  particle,  h  is  Planck's  constant  (h  =  1.055X1O34  Js),  and 

V>=4  +  4  +  A_.  (3.1.3) 

fix2  8y  dz2 


Figure  3.1  Vector  diagram  of  pairwise  electronic  charge  interactions. 
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For  electrostatic  systems,  the  potential  energy  is  generally  expressed  in  terms  of 
pairwise  interactions  between  charged  particles  (see  Fig.  3.1): 

y  =  _ L_ 

4tts0  |r2  -r, 

where  So  is  the  permittivity  of  free  space  (so  =  8.854xl0-12  C2/N-m2),  and  |  *2  - 1\  |  is  the 
distance  between  charges  q-i  and  qi. 

3.1.1  The  molecular  Hamiltonian 

Under  the  basic  assumption  of  electronic  structure  methods,  a  molecule  is  a  collection 
of  charged  quantum  particles.  The  molecular  Hamiltonian  has  the  form  of  Eq.  3.1.1 
above,  however  the  kinetic  energy  is  now  a  summation  over  all  the  particles  in  the 
molecule. 


-  l 

~  2m  ?#», 


d 2  d2  d2  ^ 


(3.1.5) 


and  the  potential  energy  is  the  Coulomb  interaction  between  each  pair  of  charged 
particles  (electron-nucleus  attraction,  nucleus-nucleus  repulsion,  and  electron-electron 
repulsion): 


v=— !— (3.1.6) 

47us0  j  k<j\rk  —  r^l 

For  electrons,  q  =  -e(e  =  1.602xlO19  C),  and  for  nuclei  of  atomic  number  Z,q  =  +Ze. 


3.2  Hartree-Fock  theory 

3.2.1  The  Schrodinger  equation 

The  quantum  mechanical  description  of  chemical  bonds  is  given  by  a  space-  and  time- 
dependent  probability  distribution:  the  molecular  wavefunction,  T'moj(t).  The  molecular 
wavefunction  is  defined  by  the  Schrodinger  equation 

H  (3.2.1) 

If  the  potential  energy  operator  is  time-independent,  then  the  solution  obtained  by 
separation  of  variables,  leads  to  the  molecular  wavefunction 
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=  (3.2.2) 

where  v| jmoi  satisfies  the  time-independent  Schrodinger  equation 

Hmo/Vmo/  =Emol\vmol,  (3.2.3) 

and  Emoi  is  the  total  energy  of  the  molecule.  Solutions  of  the  time-independent 
Schrodinger  equation  represent  various  stationary  states  of  the  molecule 
(corresponding  to  stable  or  meta-stable  electronic  configurations).  The  set  of 
wavefunctions  y  which  satisfy  Eq.  3.2.3  are  its  eigenfunctions,  and  the  energies  of  the 
molecule,  Emoi,  in  each  stationary  state  are  its  eigenvalues.  The  stationary  state  with  the 
lowest  energy  is  called  the  'ground  state/ 

3.2.2  Antisymmetry  and  electron  spin 

Standard  electronic  structure  methods  assume  that  the  molecular  wavefunction 
describing  several  electrons  can  be  written  as  a  product  of  single-electron 
wavefunctions  called  'orbitals';  that  is,  for  a  molecule  containing  n  electrons, 

W  moi  (1,2, -k)  =  M/(l)M/(2)---V|/(w) .  (3.2.4) 

Electrons  possess  an  intrinsic  angular  momentum  or  'spin'  with  a  value  of  ±%.  A  half¬ 
integer  spin  quantum  number  implies  that  electrons  are  antisymmetric  with  respect  to 
exchange  -  in  other  words,  a  wavefunction  describing  a  pair  of  electrons  i  and  j  must 
change  sign  when  the  electrons  are  interchanged: 


(3.2.5) 


The  simplest  antisymmetric  combination  of  molecular  orbitals  (MOs)  is  a  matrix 
determinant.  A  HF  wavefunction  is  constructed  by  assigning  electrons  to  molecular 
orbitals  4>(r)  in  pairs  of  opposite  spin,  and  then  forming  a  determinant  using  two 
spin  functions  a  and  p,  where 

a(T)  =  1  a(l)  =  0  (3.2.6) 

p(T)  =  0  P(40  =  l.  (3.2.7) 


For  two  electrons  i  and;,  the  total  wavefunction  takes  the  form: 


=  <Kr) 


a(i) 

at/) 


MO 

m 


(3.2.8) 


with  a  determinant 
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V(iJ)  =  ^t*(OPC/)-KOo(/)]  (3-2.9) 

which  satisfies  the  antisymmetrisation  condition  of  Eq.  3.2.4.  For  a  molecule  containing 
n  electrons,  the  wavefunction  is  referred  to  as  a  'Slater  determinant/  and  takes  the 
form: 


♦,(1)0(1)  ♦,(!»©  ♦,(1)0(1)  MW)  -  ♦^(1)0(1)  ♦j/WWl) 

♦i (2)o(2)  ♦,(2)13(2)  ♦,(2)a(2)  ♦,(2)|3(2)  -  ♦^(2)a(2)  ♦^(2)|3(2) 

♦,(«)a(«)  ♦,(«)?(«)  ♦,(n)a(«)  ♦,(n)(3(n)  —  ♦,/(")“(»)  ♦^('OPOO 

(3.2.10) 

3.2.3  Ab  initio  essentials 

For  systems  of  more  than  two  interacting  particles,  the  Schrodinger  equation  cannot  be 
solved  exactly.4  Therefore,  all  ab  initio  calculations  for  many-body  systems 
(e.g.,  molecules)  involve  some  level  of  approximation  and  indeed,  some  level  of 
empirical  parameterisation.5  Nevertheless,  ab  initio  methods  for  molecular  calculations 
must  satisfy  a  set  of  stringent  criteria: 

1.  Solutions  must  be  well  defined  and  specified  by  both  the  structure  and  the 
electronic  states  of  the  molecule. 

2.  The  potential  energy  of  the  molecule  must  vary  smoothly  and  continuously  with 
respect  to  displacements  of  the  atomic  nuclei. 

3.  The  model  must  contain  no  bias  (e.g.,  assuming  a  chemical  bond  exists  between  two 
atoms.) 

4.  The  model  must  be  'size-consistent'  -  that  is,  solutions  and  their  associated  errors 
must  scale  in  proportion  to  the  size  of  the  molecule.6 

5.  The  model  must  be  'variational'  -  that  is,  approximate  solutions  must  provide  an 
upper  bound  to  the  'true'  energy  of  the  system.  Consequently,  the  approximate 
solution  having  the  lowest  energy  represents  the  closest  fit  to  the  true 
wavefunction,  within  the  constraints  of  the  method. 


4  This  is  true  for  any  system,  quantum  or  classical,  since  current  permutation  operators  (e.g.,  the 
Levi-Civita  symbol,  sijk)  can  only  account  for  pairwise  (two-particle)  relationships  [13]. 

5  Indeed,  ab  initio  electronic  structure  calculations  could  be  considered  the  'ultimate'  semi- 
empirical  approximation!  It  is  perhaps  more  consistent  to  describe  ab  initio  methods  as 
'universal,'  that  is,  the  a  set  of  wavefunctions  are  optimised  for  a  particular  method  using 
empirical  data  from  a  set  of  atoms  and  molecules,  and  it  is  then  assumed  that  these 
wavefunctions  may  be  applied  to  calculations  involving  all  atoms  and  molecules. 

6  Size-consistency  allows  direct  comparison  between  molecules  of  different  sizes,  and  is 
therefore  an  important  consideration  for  studies  of  chemical  reactions. 
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3.2.4  Approximations 

3. 2.4.1  Non-relativistic  solutions 


Typical  ab  initio  electronic  structure  methods  are  time-independent,  solving  only  the 
spatial  wavefunction,  iy(r).  This  approximation  is  reasonable  for  most  chemical 
systems,  but  breaks  down  in  calculations  involving  large  atoms  (late  transition  metals, 
lanthanides,  etc.),  whose  inner-shell  electrons  have  velocities  approaching  the  speed  of 
light 

32.4.2  Bom-Oppenheimer  approximation 

Electrons  in  molecules  are  much  lighter  than  nuclei,  and  therefore  generally  have  much 
higher  velocities.  Hence,  under  most  circumstances,  one  can  assume  that  electrons 
respond  instantaneously  to  nuclear  displacements.  In  practice,  this  means  that  the 
molecular  Hamiltonian  can  be  written  assuming  the  nuclear  positions  are  fixed  (i.e., 
neglecting  nuclear  kinetic  energy  terms): 


H 


mol 


:  T  +V 

1  elec  ^  y 


_ ij 2  electrons 

~z~~  I  V?  + 

i 


1  ^  V  Qflk 

47TS0  j  k>j  rk  -  Tj 


(3.2.11) 


Most  ab  initio  calculations  solve  only  the  electronic  part  of  the  molecular  wavefunction, 
and  therefore  cannot  account  for  systems  where  the  electronic  states  are  strongly 
coupled  to  nuclear  vibrations. 


32.4.3  Single  particle  approximation 

As  mentioned  in  Section  3.2.2,  standard  electronic  structure  methods  approximate  the 
total  wavefunction  of  a  many-electron  system  as  the  product  of  single-electron 
wavefunctions  (Eq.  3.2.4).  This  is  the  essence  of  Hartree-Fock  (HF)  theory,  which 
describes  each  electron  in  a  molecule  as  moving  in  the  average  electric  field  generated 
by  the  other  electrons  and  nuclei.  As  a  single-particle  theory,  HF  theory  systematically 
overestimates  molecular  energies  because  it  neglects  the  correlated  motion  of  electrons 
resulting  from  Coulomb  interactions. 

32.4.4  Linear  combinations  of  atomic  orbitals  (LCAO) 

Although  there  is  no  exact  analytical  solution  to  the  time-independent  molecular 
SchrSdinger  equation  (Eq.  3.2.3)  for  systems  containing  more  than  one  electron, 
approximate  solutions  can  be  obtained  using  standard  numerical  techniques.  The 
approach  of  all  ab  initio  techniques  is  to  build  the  total  wavefunction  from  a  'basis'  set 
of  mathematical  functions  capable  of  reproducing  critical  properties  of  the  system.  An 
individual  molecular  orbital  may  then  be  expressed  as 
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i(r)  =  ZcA(r),  (3.2.12) 

where  ^(r)  are  the  basis  functions,  and  the  coefficients  c„;  are  adjustable  parameters. 
For  a  molecular  wavefunction,  the  electronic  orbitals  of  the  constituent  atoms  form  a 
natural  set  of  basis  functions.  These  atomic  orbitals  can  in  turn  be  represented  by 
different  types  of  mathematical  functions.  A  highly  accurate  set  of  atomic  orbitals 
(Slater-type  orbitals  or  STOs)  are  based  on  hydrogenic  wavefunctions  having  the  form 

Xsro(r)  ~  Ce~ar .  (3.2.13) 

Exponential  functions  are  not  well  suited  to  numerical  manipulation,  so  most  electronic 
structure  calculations  approximate  STOs  with  a  linear  combination  of  gaussian-type 
functions. 


Xsro(r)  X(i  —  ' 

V 


(3.2.14) 


where  and  av  are  adjustable  parameters.  As  can  be  seen  from  Fig.  3.2,  gaussian-type 
functions  provide  reasonable  approximations  of  STOs,  except  at  very  small  or  very 
large  electron-nucleus  separations. 


Figure  3.2  (a)  Comparison  of  exponential  and  gaussian  functions,  (b)  Comparison  of  the  same 
exponential  function  and  a  sum  of  three  gaussians. 

Linear  combinations  of  'primitive'  gaussian  functions  are  referred  to  as  'contracted' 
gaussians.  Standard  ah  initio  software  packages  like  Gaussian98  [50]  offer  a  choice  of 
basis  sets  containing  contracted  gaussians  optimised  to  reproduce  the  chemistry  of  a 
large  range  of  molecular  systems.  A  more  detailed  discussion  of  molecular  basis  sets  is 
provided  in  section  3.3. 


13 


DSTO-GD-0253 


3.2.5  The  variational  principle  and  Roothan-Hall  equations. 

The  sections  above  describe  a  method  for  constructing  a  determinantal  wavefunction 
from  a  set  of  LCAO-MOs.  It  remains  now  to  define  a  method  to  determine  the  MO 
expansion  coefficients  c«  (Eq.  3.2.11)  which  optimise  the  molecular  wavefunction  for  a 
particular  system.  The  variational  nature  of  the  model  (Section  3.2.3)  guarantees  that 
the  energy  eigenvalue  solution  for  any  approximate  wavefunction  is  always  greater 
than  the  energy  obtained  using  the  exact  wavefunction.  It  follows  that  the  set  of 
coefficients  which  minimise  the  energy  of  the  resultant  wavefunction  will  give  the  best 
approximation  to  the  exact  wavefunction  from  a  chosen  basis  set. 

The  variational  constraint  leads  to  a  set  of  algebraic  equations  (Roothan-Hall)  for  cw, 
expressed  in  matrix  form  as 


FC  =  SCs 


(3.2.15) 


where 

•  C  is  die  matrix  of  MO  expansion  coefficients; 

•  F  is  the  Fock  matrix,  which  is  the  sum  of  a  term  representing  the  energy  of  a 
single  electron  in  the  field  of  the  bare  atomic  nuclei  and  a  term  describing 
electron-electron  repulsion  within  an  averaged  field  of  electron  density; 

•  S  is  a  matrix  describing  the  overlap  of  molecular  orbitals;  and 

•  e  is  a  diagonal  matrix  containing  the  one-electron  energies  of  each  molecular 
orbital 

Since  the  terms  within  the  Fock  matrix  F  depend  upon  the  electron  density,  which  in 
turn,  depends  upon  molecular  wavefunction  defined  by  the  matrix  of  MO  expansion 
coefficients  C,  the  Roothan-Hall  equations  are  nonlinear,  and  must  be  solved  by  an 
iterative  procedure  termed  the  'self-consistent  field'  (SCF)  method.  Upon  convergence 
of  the  SCF  method,  the  minimum-energy  MOs  produce  the  electric  field  which 
generate  the  same  orbitals  (hence,  the  self-consistency). 


3.3  Basis  Sets 

In  general,  a  basis  set  is  an  assortment  of  mathematical  functions  used  to  solve  a 
differential  equation.  In  quantum  chemical  calculations,  the  term  'basis  set'  is  applied 
to  a  collection  of  contracted  gaussians  representing  atomic  orbitals,  which  are 
optimised  to  reproduce  the  desired  chemical  properties  of  a  system. 

Standard  ab  initio  software  packages  generally  provide  a  choice  of  basis  sets  that  vary 
both  in  size  and  in  their  description  of  the  electrons  in  different  atomic  orbitals.  Larger 
basis  sets  include  more  and  a  greater  range  of  basis  functions.  Therefore,  larger  basis 
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sets  can  better  refine  the  approximation  to  the  'true'  molecular  wavefunction7,  but 
require  correspondingly  more  computer  resources.  Alternatively,  accurate 
wavefunctions  may  be  obtained  from  different  treatments  of  electrons  in  atoms.  For 
instance,  molecules  containing  large  atoms  (Z  >  30)  are  often  modelled  using  basis  sets 
incorporating  approximate  treatments  of  inner-shell  electrons  which  account  for 
relativistic  phenomena. 

'Minimal'  basis  sets  contain  the  minimum  number  of  AO  basis  functions  needed  to 
describe  each  atom  ( e.g .,  Is  for  H  and  He;  Is,  2s,  2px,  2py,  2pz  for  Li  to  Ne).  An  example 
of  a  minimal  basis  set  is  STO-3G,  which  uses  three  gaussian-type  functions  (3G)  per 
basis  function  to  approximate  the  atomic  Slater-type  orbitals  {see  Fig.  3.2b).  Although 
minimal  basis  sets  are  not  recommended  for  consistent  and  accurate  predictions  of 
molecular  energies,  their  simple  structure  provides  a  good  tool  for  visualising 
qualitative  aspects  of  chemical  bonding.  Improvements  on  minimal  basis  sets  are 
described  below  and  illustrated  in  Fig.  3.3. 

3.3.1  Split  valence  basis  sets 

In  split  valence  basis  sets,  additional  basis  functions  (one  contracted  gaussian  plus 
some  primitive  gaussians)  are  allocated  to  each  valence  atomic  orbital.  The  resultant 
linear  combination  allows  the  atomic  orbitals  to  adjust  independently  for  a  given 
molecular  environment.  Split  valence  basis  sets  are  characterised  by  the  number  of 
functions  assigned  to  valence  orbitals.  'Double  zeta'  basis  sets  use  two  basis  functions 
to  describe  valence  electrons,  'triple  zeta'  use  three  functions,  and  so  forth.  Basis  sets 
developed  by  Pople  and  coworkers  [51]  are  denoted  by  the  number  of  gaussian 
functions  used  to  describe  inner  and  outer  shell  electrons.  Thus  '6-21G'  describes  an 
inner  shell  atomic  orbital  with  a  contracted  gaussian  composed  of  six  primitive 
gaussians,  an  inner  valence  shell  with  a  contracted  gaussian  composed  of  two 
primitives,  and  an  outer  valence  shell  with  one  primitive.  Other  split- valence  sets 
include  3-21G,  4-31G,  and  6-311G. 


Figure  3.3.  Basis  set  improvements. 


7  It  must  be  emphasised  that  although  MOs  are  used  to  obtain  values  for  physical  'observables' 
(e.g.,  molecular  heats  of  formation,  electronic  dipole  moments),  they  are  strictly  mathematical 
constructs  (indeed,  they  are  solutions  to  the  Hartree-Fock  approximation  as  opposed  to  the  true 
many-body  electronic  Hamiltonian),  and  do  not  themselves  contain  any  physical  meaning. 
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3.3.2  Polarised  basis  sets 

Polarisation  functions  can  be  added  to  basis  sets  to  allow  for  non-uniform 
displacement  of  charge  away  from  atomic  nuclei,  thereby  improving  descriptions  of 
chemical  bonding.  Polarisation  functions  describe  orbitals  of  higher  angular 
momentum  quantum  number  than  those  required  for  the  isolated  atom  (e.g.,  p-type 
functions  for  H  and  He,  and  d-type  functions  for  atoms  with  Z  >  2),  and  are  added  to 
the  valence  electron  shells.  For  example,  the  6-31G(d)  basis  set  is  constructed  by  adding 
six  d-type  gaussian  primitives  to  the  6-31G  description  of  each  non-hydrogen  atom. 
The  6-31G(d,p)  is  identical  to  6-31G(d)  for  heavy  atoms,  but  adds  a  set  of  gaussian  p- 
type  functions  to  hydrogen  and  helium  atoms.  The  addition  of  p-orbitals  to  hydrogen  is 
particularly  important  in  systems  where  hydrogen  is  a  bridging  atom. 

3.3.3  Diffuse  basis  sets 

Species  with  significant  electron  density  far  removed  from  the  nuclear  centres 
(e.g.,  anions,  lone  pairs  and  excited  states)  require  diffuse  functions  to  account  for  the 
outermost  weakly  bound  electrons.  Diffuse  basis  sets  are  recommended  for 
calculations  of  electron  affinities,  proton  affinities,  inversion  barriers  and  bond  angles 
in  anions.  The  addition  of  diffuse  s-  and  p- type  gaussian  functions  to  non-hydrogen 
atoms  is  denoted  by  a  plus  sign  -  as  in  '3-21+G/  Further  addition  of  diffuse  functions 
to  both  hydrogen  and  larger  atoms  is  indicated  by  a  double  plus8. 

3.3.4  High  angular  momentum  basis  sets 

Basis  sets  with  multiple  polarisation  functions  are  now  practical  for  many  systems  and, 
although  not  generally  required  for  Hartree-Fock  calculations,  are  useful  for  describing 
the  interactions  between  electrons  in  electron  correlation  methods  (see  section  3.4). 
Examples  of  high  angular  momentum  basis  sets  include: 

•  6-31G(2d)  -  two  d-functions  are  added  to  heavy  atoms; 

•  6-311G(2df,  pd)  -  besides  the  (311)  valence  functions,  two  d  functions  and  one  / 
function  are  added  to  heavy  atoms,  and  p  and  d  functions  to  hydrogen; 

•  6-311G(3df,  2df,  p)  -  three  d  functions  and  one  f  function  are  added  to  atoms  with 
Z  >  11,  two  d  functions  and  one /function  to  first-row  atoms  (Li  to  Ne)  and  one  p 
function  to  hydrogens. 

High  angular  momentum  basis  sets  augmented  with  diffuse  functions  represent  the 
most  sophisticated  basis  sets  available.  As  is  discussed  in  Section  4.3.3  below,  energetic 
molecules  have  a  complicated  electronic  structure  due  to  the  presence  of  nitro  groups. 
Therefore,  the  most  accurate  ab  initio  studies  of  energetic  materials  would  be  produced 
by  reasonably  sophisticated  polarised  split-valence  basis  sets  augmented  with  high 
angular  momentum  and  diffuse  atomic  orbitals.  However,  the  size  of  the  optimum 
basis  set,  especially  when  used  with  electron  correlation  methods  (discussed  in 


8  Diffuse  functions  on  hydrogen  atoms  seldom  make  a  difference  in  accuracy  [1], 
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Section  3.3  below)  will  ultimately  be  determined  by  the  size  of  the  energetic  molecule, 
the  amount  computing  power  available,  and  the  time  allotted  for  the  studies.  For 
example,  conformational  studies  of  RDX  (1,3,5  trinitro-l,3^-triazacyclohexane)  range 
from  rudimentary  HF/STO-3G  calculations  performed  in  1991  [52]  to  more 
sophisticated  investigations  in  1997  using  B3LYP  density  functional  theory  and 
6-311G(d,p)  [53]. 


3.4  Electron  correlation 

Hartree-Fock  theory  is  a  single-particle  approximation,  and  therefore  cannot 
adequately  treat  the  correlated  motion  of  electrons  that  occurs  due  to  electron-electron 
interactions.  Neglect  of  electron  correlation  has  been  blamed  for  systematic  FIF  errors 
such  as  underestimated  bond  lengths  and  overestimated  vibrational  frequencies. 
Calculations  added  to  HF-SCF  theory  to  remedy  these  errors  are  termed  'electron 
correlation'  or  'post-HF'  methods.  There  are  three  general  types  of  electron  correlation 
treatments:  configuration-interaction  (Cl)  methods,  Moller-Plesset  (MP)  perturbation 
theory,  and  density  functional  theory  (DFT). 

3.4.1  Configuration-interaction 

A  HF  wavefunction  of  a  molecule  has  one  determinant;  therefore,  it  can  only  describe 
only  one  electronic  configuration  in  a  molecule.  This  approximation  is  so  restrictive 
that  HF  theory  cannot  predict  dissociation  of  the  simplest  molecule,  the  hydrogen 
dimer,  into  its  neutral  atomic  components.  However,  correct  dissociation  behaviour  is 
restored  when  extra  configurations,  corresponding  to  electronically  excited  states,  are 
added  to  the  wavefunction.  Configuration-interaction  methods  incorporate  excited 
state  configurations  into  the  wavefunction  by  constructing  new  determinants  from  the 
original  HF  determinant.  New  determinants  are  created  by  replacing  one  or  more 
occupied  orbitals  with  unoccupied  (virtual)  orbitals  of  higher  energy.  The  number  of 
replacements  within  the  determinants  designates  the  level  of  Cl.  For  instance,  single 
substitution  (CIS)  switches  one  pair  of  occupied  and  virtual  orbitals, 

HF  ~  l^l, ---4)  j’-§HOMO’§LUMO->—§i’—§i\ 

Vc/S  =|<t)l,— 1 HOMO’§  LUMO>"-§  j  ’"•<t)n|  (3.4.1) 

and  is  equivalent  to  a  one-electron  excitation.  Higher-order  calculations  include  CID 
(double  substitutions),  which  generates  determinants  where  two  orbital  pairs  are 
switched;  CISD,  which  adds  singly  and  double-substituted  determinants;  and  CISD(T), 
with  single,  double  and  triple  excitations.  The  theoretical  limit  of  this  expansion  -  a  full 
Cl  calculation  -  forms  the  molecular  wavefunction  as  a  linear  combination  of  the  HF 
determinant  and  all  possible  substituted  determinants: 
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V fiiii-ci  ~  ao Vhf  +  Sfln¥s,z),r,...  (3.4.2) 

n>  0 

where  the  coefficients  a,  are  determined  by  minimising  the  energy  of  the  total 
wavefunction.  Full  Cl  provides  the  most  complete  non-relativistic  treatment  possible 
for  a  molecular  system;  however,  it  is  extremely  computer-intensive,  and  generally 
untenable  for  all  but  the  simplest  molecules  described  by  small  basis  sets. 

3.4.2  M0ller-Plesset  perturbation 

For  typical  molecules,  most  of  the  ground-state  energy  originates  from  one-electron  HF 
contributions.  Moller-Plesset  perturbation  theory  assumes  that  the  effects  of  electron 
correlation  are  minor,  and  can  be  described  by  small  corrections  (perturbations)  to  the 
HF  solution.  Essentially,  MP  methods  assume  the  true  molecular  Hamiltonian  (3.2.3) 
can  be  divided  into  two  parts: 


Hm0/  -  Hone_e+ A,Pman>,_e  (3.4.3) 

where  Hone_e  denotes  single-electron  energy  contributions  that  can  be  solved  exactly  by 

HF-SCF,  and  Pmflny_e  represents  contributions  due  to  electron  correlation.  The 

coefficient  X  is  used  to  generate  power  series  expansions  of  the  energy  and  the 
molecular  wavefunction: 

Emol  =  E(0)  +  XEm  +  X2E(2)  +  X3E{3)  + ...  (3.4.4) 

\ymol  =  i{/(0)  +  ^y(1)  +  ^.y 2)  +  y3)  + ...  (3.4.5) 

These  series  are  then  substituted  back  into  the  molecular  Schr 6 dinger  equation  with 
the  modified  Hamiltonian  (3.4.3),  and  each  term  is  evaluated  in  turn.  The  first  two 
energy  terms  constitute  the  HF  energy 

EHF=Em+Em,  (3.4.6) 

and  are  used  to  evaluate  vj/t1),  which  is  composed  of  a  linear  combination  of  substituted 
determinants 


V')=Zbw,'  (3A7) 

n 

with  coefficients  bi  that  are  inversely  proportional  to  the  difference  in  energy  between 
the  ground  state  and  the  associated  excited  states.  That  is,  substituted  wavefunctions 
close  in  energy  to  the  ground  state  make  larger  contributions  to  the  perturbation 
expansion.  The  first-order  correction  to  the  wavefunction,  if/.1),  is  then  used  to  calculate 
the  second-order  correction  to  the  total  energy  E(2),  and  so  on,  back  and  forth,  until  the 
desired  order  of  correction  is  achieved.  Moller-Plesset  theory  is  designated  by  the  order 


18 


DSTO-GD-0253 


of  the  perturbative  corrections.  Since  the  second-order  energy  term  E(2)  is  the  first-order 
correction  to  the  HF  energy,  the  designation  begins  at  MP2,  and  continues  with  MP3, 
MP4,  MP5,  and  so  forth. 

3.4.3  Density-functional  theory 

DFT  theory  models  electron  correlation  as  a  functional9  of  the  electron  density,  p.  The 
functional  employed  by  current  DFT  methods  partitions  the  electronic  energy  via  the 
Kohn-Sham  equations  [54,55]  into  several  terms 

E  =  ET +EV +EJ +EXC  (3.4.8) 

where  ET  is  the  kinetic  energy  term  (arising  from  the  motion  of  the  electrons),  Ev  is  the 
potential  energy  term  that  includes  nuclear-electron  and  nuclear-nuclear  interactions,  B 
is  the  electron-electron  repulsion  term  and  Exc  is  the  electron  correlation  term.  All  terms 
except  nuclear-nuclear  repulsions  are  functions  of  the  electron  density.  The  terms  ET  + 
Ev  +  B  represent  the  classical  energy  of  the  electron  distribution,  while  Exc  represents 
both  the  quantum  mechanical  exchange  energy,  which  accounts  for  electron  spin10,  and 
the  dynamic  correlation  energy  due  to  the  concerted  motion  of  individual  electrons. 
Pure  DFT  methods  calculate  Exc  by  pairing  an  exchange  functional  with  a  correlation 
functional  and  so  are  designated  by  the  choice  of  combination.  For  example,  BLYP 
combines  Becke's  gradient-corrected  exchange  functional  with  the  gradient-corrected 
correlation  functional  of  Lee,  Yang  and  Parr  [56]. 

DFT  calculations  fall  into  three  general  categories:  local  density  approximations  (LDA), 
generalised  gradient  approximations  (GGA),  and  'hybrid'  combinations  of  DFT  and 
Hartree-Fock  terms.  LDA  exchange  and  correlation  functionals  only  contain  terms 
related  to  electron  density11  -  an  approach  that  works  for  some  bulk  materials,  but  fails 
to  accurately  predict  properties  in  isolated  molecules.  GGA  ('nonlocal')  functionals 
contain  terms  that  depend  upon  both  die  electron  density  and  the  density  gradients12. 


9  Afunctional  is  a  function  defined  by  another  function,  that  is,  a  function  of  a  function. 

10  The  quantum-mechanical  exchange  energy  is  accounted  for  in  HF  theory  by  generating 
antisymmetric  determinental  wavefuntions,  as  discussed  in  section  2.2.2.  Thus,  Coulomb  "self¬ 
interaction"  energy  can  be  considered  a  type  of  electron  correlation  already  included  in  HF 
theory. 

11  An  example  of  a  local  functional  is  the  'Slater'  functional: 

V 

where  p,  the  electron  density,  is  a  function  of  r .  This  local  functional  reproduces  the  properties 
of  a  uniform  electron  gas.  Note:  this  is  not  the  same  as  Slater  determinants  discussed  in  the 
previous  footnote. 

12Gradients  represent  spatial  changes  of  a  function  as  can  be  quite  large.  For  example,  the 
electronic  density  within  a  metal  might  be  quite  uniform,  and  thus  be  well  modelled  by  the 
Slater  functional  in  footnote  6,  yet  display  large  gradients  at  the  metal  surface.  Large  gradients 
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The  gradient-corrected  density  functional  method  BLYP  is  capable  of  predicting 
intramolecular  bond  dissociation  energies  to  within  a  few  kj/mol  [57].  However,  the 
generalised  gradient  approximation  severely  underestimates  activation  barriers  for 
some  reactions  due  to  neglect  of  Coulomb  'self-interaction7  of  the  electrons  [58].  This 
problem  is  remedied  with  hybrid  methods  that  combine  Hartree-Fock  self-interaction 
corrections  with  density  functional  exchange  and  correlation.  Examples  of  hybrid 
methods  are  B3LYP  and  B3PW91,  where  B3  denotes  Becke's  three-parameter  hybrid 
functional  [59],  while  'PW9T  and  'LYF  are  gradient-corrected  correlation  functionals  of 
Perdew  and  Wang  [60]  and,  as  above,  Lee,  Yang  and  Parr. 

3.4.4  Caveats  for  post-HF  methods 

Although  post-SCF  treatments  should  improve  upon  Hartree-Fock  theory,  each  of 
these  techniques  violates  at  least  one  of  the  fundamental  criteria  that  define  ab  initio 
methodology.  Truncated  Cl  calculations  (i.e.,  anything  less  than  full  configuration 
interaction)  are  not  size-consistent,  so  comparative  Cl  studies  should  include  added 
corrections13.  MP  perturbation  theory  is  not  a  variational  technique,  and  as  such  does 
not  provide  an  upper  bound  to  the  true  molecular  energy14.  Finally,  DFT  contains  a 
bias:  the  exact  functional  dependence  of  Exc  upon  the  electron  density  is  unknown  and 
must  be  approximated  by  assuming  that  p  has  certain  properties  (e.g.,  it  behaves  like  a 
uniform  electron  gas).  Excited-state  DFT  calculations,  which  require  the  application  of 
linear  response  theory,  has  so  far  been  limited  to  atomic  systems15.  These  caveats  serve 
as  a  reminder  that  although  ab  initio  methods  are  capable  of  yielding  predictions  close 
to  experimental  results,  many  of  the  approximations  that  make  quantum  chemical 
computations  feasible  are  wholly  nonphysical.  Hence,  even  the  most  sophisticated 
ab  initio  software  packages  should  not  be  treated  simply  as  'black  boxes.' 

Given  these  considerations,  the  results  of  any  MO  calculation  on  energetic  materials  are 
subject  to  scrutiny.  Nevertheless,  one  can  assume  with  a  high  level  of  confidence  that 
DFT/6-31G(d)  will  yield  bond  lengths  accurate  to  +0.01  A,  and  the  molecular 
vibrational  frequencies  to  within  +100  cm-1.  However,  DFT  predictions  of  chemical 
reaction  energies  should  be  treated  as  qualitatively  correct  at  best.  For  quantitative 


also  exist  across  interfaces  of  p-  and  n-  type  semiconductors,  and  govern  gating  and  switching 
of  electric  current  in  transistors. 

13  One  method  ,  quadratic  configuration  interaction  (QCI),  restores  size  consistency  by  adding 
terms,  with  corrections  available  up  to  quadruply  excited  states:  QCISD,  QCISD(T)  and 
QCISD(TQ).  Coupled  cluster  methods  correct  for  size-consistency  by  accounting  for 
simultaneous  substitutions  in  different  molecules  -  a  technique  useful  for  modelling  chemical 
reactions.  Coupled-cluster  calculations  (CCD,  CCSD)  are  useful  when  electron  correlation 
makes  a  major  contribution  to  electron-pair  bonds  and  so  cannot  be  reasonably  described  by 
perturbation  methods  [61]. 

14  The  second-order  term  (MP2)  is  always  negative  and  often  overcorrects  for  electron 
correlation. 

15  In  other  words,  DFT  should  not  be  used  to  calculate  excited  states  of  molecules  with  existing 
commercial  software  [53], 
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accuracy  (±  40  kj / mol),  a  recommended  procedure  is  to  optimise  the  molecular 
structure  using  DFT  with  a  high-level  split-valence  polarised  basis  set,  then  recalculate 
the  structural  energy  with  MP4  and  the  same  or  a  slightly  less  sophisticated  split- 
valence  polarised  basis  set. 


3.5  Open-shell  calculations  and  spin  contamination 


In  addition  to  the  considerations  above,  input  to  molecular  orbital  calculations  must 
specify  the  treatment  of  electron  spin  -  that  is,  whether  the  molecule  will  be  described 
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Figure  3.4.  Electron  assignments  for  closed-shell  and  open-shell  calculations. 

using  an  'open-shell'  or  'closed-shell'  model.  Most  molecules  have  singlet  electronic 
ground  states  (zero  total  electron  spin):  they  contain  an  even  number  of  electrons 
which  can  be  assigned  to  orbitals  in  pairs  of  opposite  spin.  HF  models  that  pair 
electrons  in  this  way  are  referred  to  as  'spin-restricted  closed-shell'  calculations,  and 
are  adequate  for  treating  most  types  of  chemical  systems.  Flowever,  there  are  several 
cases  which  require  explicit  consideration  of  unpaired  electrons,  and  are  better 
described  by  open-shell  unrestricted  models  that  generate  separate  MOs  for  spin-up 
(a)  and  spin-down  ((3)  electrons.  Spin-unrestricted  calculations  are  used  for: 


•  Molecules  with  odd  numbers  of  electrons  ( e.g .,  some  ions  and  radicals); 

•  Excited  states; 

•  Molecules  with  triplet  or  higher-spin  ground  states  (e.g.,  O2); 

•  Processes  that  separate  electron  pairs  (e.g.,  dissociation);  and 

•  Delocalised  orbitals  for  resonant  systems. 


Spin-unrestricted  wavefunctions  are  not  eigenfunctions  of  the  spin  operator  (S  ),  and 
therefore  suffer  from  intermixing  between  spin  states.  In  most  open-shell  systems,  spin 
contamination  is  negligible;  however,  for  molecules  where  the  unpaired  electrons  are 
delocalised,  contamination  can  be  quite  significant,  especially  from  the  next-higher 
spin  state.  Unrestricted  calculations  of  nitric  oxide,  a  diatomic  molecule  with  one 
impaired  electron  and  a  doublet  (2n)  ground  state,  display  extreme  spin  contamination 
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when  small  basis  sets  (STO-3G  or  3-21G)  are  used16.  Spin  contamination  can  also  be 
expected  in  calculations  of  energetic  materials  containing  nitro  (-NO2)  groups 
involving  excited  states  or  radical  producing  reactions. 

The  detrimental  effects  of  spin  contamination  can  be  appreciable,  although  the  degree 
of  spin  contamination  decreases  with  increasing  basis  set  size.  Spin  contamination  has 
been  blamed  for  poor  predictions  of  molecular  geometries  [62]  and  vibrational 
frequencies  [63].  The  effects  of  spin  contamination  are  particularly  acute  when 
occupied  and  unoccupied  orbitals  are  nearly  degenerate.  This  can  happen,  for  instance, 
when  a  bond  is  stretched  to  the  limit  of  dissociation,  and  the  separation  between 
bonding  and  antibonding  orbitals  approaches  zero,  generating  a  wavefunction  that 
contains  an  equal  mixture  of  singlet  and  triplet  states.  Another  example  is  ethylene 
torsion,  where  the  n  and  71*  orbitals  become  equivalent  by  symmetry  at  the  transition 
state  for  rotation  about  the  double  bond  [64]. 

It  is  more  desirable  to  remove  spin  contamination  than  to  ignore  it,  however,  spin 
annihilation  is  beyond  the  scope  of  simpler  post-HF  methods.  Open-shell  HF 
wavefunctions  effectively  introduce  some  electron  correlation,  because  the  use  of 
different  orbitals  for  each  electron  allows  a  and  (3  electrons  to  be  spatially  further  apart. 
Hence,  simpler  correlation  corrections  (MP2  and  CIS)  do  not  generally  improve  spin- 
contaminated  HF  wavefunctions  [65].  Suggested  spin  annihilation  methods  include 
high  levels  of  electron  correlation,  spin-projected  energy  corrections,  or  the  use  of 
various  forms  of  restricted  open-shell  wavefunctions.  Gaussian  98  automatically 
incorporates  annihilation  of  one  spin  contaminant  into  its  UHF  calculations.  Density 
functional  methods  also  appear  to  reduce  spin  contamination.  Baker,  Scheiner  and 
Andzelm  [66]  observed  that  'pure'  unrestricted  Kohn-Sham  wavefunctions  display 
significantly  less  spin  contamination  than  their  HF  counterparts,  and  concluded  that 
spin  contamination  arises  primarily  from  the  use  of  Slater  determinants  in  HF  theory. 
Therefore,  if  spin  contamination  is  a  consideration  in  density  functional  calculations, 
hybrid  methods  that  include  HF  exchange  should  be  avoided. 


16  An  obvious  manifestation  of  spin  contamination  is  an  overestimate  of  spin-squared 
expectation  values.  For  a  spin  eigenfunction  4>(a,P),  S2  <{>(a,P)  =  s(s+l)4>(a,p),  where  s  is  the  total 
spin.  For  pure  doublet  states  (s  =  V2),  <S2>  -V2  (V2  +  1)  =  0.75.  As-mentioned  unrestricted  HF 
calculations  of  nitric  oxide  yield  values  of  <  S2  >  ranging  from  1.0  to  1.45. 
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4.  Accuracy  versus  Efficiency  for  Ab  initio  Methods 

4.1  Error  reduction 

4.1.1  Methods 

The  standard  MO  treatment  for  most  closed-shell  molecules  is  a  spin-restricted 
Hartree-Fock  self-consistent  field  calculation.  HF-SCF  calculations  generally  yield 
accurate  molecular  structures  but  are  less  successful  at  predicting  molecular  energies. 
Depending  upon  the  system  and  the  choice  of  basis  set,  HF  energies  typically  fall 
anywhere  within  200  kj/mol  of  experimental  values,  although  errors  can  be  as  large  as 
700  kj/mol.  The  main  source  of  error  in  HF  calculations  is  neglect  of  electron 
correlation,  which  results  in  systematic  overestimates  of  molecular  energies.  To 
account  for  electron  correlation,  calculations  are  added  to  HF  results,  including 
perturbation  treatments  (e.g.,  MP2  and  MP4),  multiconfiguration  methods  (e.g.,  CISD, 
CCSD  and  CASSCF17),  and  density  functional  theory  (e.g.,  BLYP,  and  B3PW91). 
Perturbation  methods  and  density  functional  theory  (DFT)  are  typically  accurate  to 
within  40  kj/mol,  while  the  results  of  multiconfiguration  treatments  are  even  better, 
falling  within  5-8  kj/mol  of  experiment.  Since  multiconfiguration  calculations  are 
extremely  expensive  for  even  moderately  sized  molecules,  perturbation  treatments  and 
DFT  tend  to  be  the  methods  of  choice  for  treating  electron  correlation. 

4.1.2  Basis  sets 

In  general,  ab  initio  MO  studies  on  complex  unknown  systems  should  begin  with 
calculations  using  small  basis  sets,  both  to  obtain  a  timely  qualitative  assessment  of 
molecular  properties,  and  to  determine  the  computational  cost  of  the  higher-order 
calculations.  However,  minimal  basis  sets  like  STO-3G  should  only  be  used  when  the 
size  of  the  molecule  is  beyond  the  scope  of  available  computer  resources  for  even  small 
split-valence  sets  like  3-21G.  Quantitative  accuracy  improves  with  the  size  of  basis  sets, 
since  larger  basis  sets  contain  more  adjustable  parameters  and  are  thus  better  able  to 
approximate  true  molecular  wavefunctions.  Polarisation  functions  improve  the 
accuracy  of  fundamental  basis  sets  (e.g.,  changing  6-31G  to  6-31G(d)),  although  the 
ultimate  gain  is  limited18.  Diffuse  functions  (e.g.,  changing  6-31G  to  6-31+G)  will 
improve  results  under  certain  circumstances.  For  general  use,  the  smallest  standard 
basis  set  recommended  by  developers  of  Gaussian  98  is  6-31G(d)  [51]. 


17  CASSCF  is  the  acronym  for  'complete  active-space  self-consistent  field'  calculations.  CASSCF 
is  a  'full'  configuration-interaction  calculation  over  a  limited  number  of  electrons.  That  is.  Cl 
determinants  are  constructed  for  specific  molecular  orbitals,  usually  those  containing  valence 
electrons  [51]. 

18  That  is,  for  a  given  property,  calculations  performed  with  the  6-311+G(d,p)  may  be  just  as 
accurate  as  those  using  6-311++G(3df,3pd),  in  which  case  6-311+G(d,p)  defines  the  'basis  set 
limit.' 
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4.2  Resource  usage 

Given  the  limits  of  computer  memory  and  speed,  the  size  of  the  chemical  system 
generally  determines  the  maximum  size  of  the  basis  set  used  with  MO  methods.  For 
instance,  a  moderately  sized  energetic  molecule,  1,3,3-trinitroazetidine  (TNAZ), 
contains  49  pairs  of  electrons.  The  6-31G(d)  basis  set  assigns  203  basis  functions  to  the 
molecule,  49  of  which  are  occupied,  and  the  remaining  154,  unoccupied  or  'virtual' 
orbitals.  To  estimate  memory  requirements  (in  8-byte  words),  Gaussian  98  references 
[51]  suggest  using 

M  +  2Nb2  (4.1.1) 

where  Nb  is  the  number  of  basis  sets  used  in  the  calculation  and  the  value  of  M 
depends  upon  the  method  of  calculation  (typically  around  4-6  MW). 

Table  4.1  Resource  usage  by  method  and  basis  set 


Method 

CPU 

Memory 

Disk 

SCF 

N4 

N2 

N4 

MP2 

ON4 

N2 

N4 

MP4/QCISD(T) 

Q3V4 

N2 

N4 

N  =  total  number  of  basis  functions 
O  =  total  number  of  occupied  orbitals 
V  =  total  number  of  virtual  orbitals 


According  to  Eq.  4.1.1,  a  Hartree-Fock  SCF  single-point  energy  calculation  using  the 
6-31G(d)  basis  set  (denoted  'FIF/6-31G(d)')  requires  approximately  4.3  MW  (34  MB) 
worth  of  memory.  On  a  SGI  Origin  200019,  the  recorded  CPU  time  for  this  calculation  is 
19  minutes.  From  the  guidelines  listed  in  Table  4.1,  a  calculation  using  a  larger  basis  set 
will  require  significantly  more  memory  and  CPU  time.  For  instance,  a  HF/aug-cc- 
pVTZ20  SPE,  which  assigns  609  basis  functions  to  TNAZ,  would  require  at  least  40MB 
of  memory,  and  take  over  a  day  to  run21.  Electron  correlation  treatments  are  also 
expensive:  a  MP2/6-31G(d)  SPE  would  take  roughly  16  hours  to  run,  while  MP4  or 
QCISD(T)/6-31G(d)  would  take  a  year  and  a  half.  Clearly,  sophisticated  post-HF 
calculations,  such  as  Cl  and  MP2  coupled  with  correlation-consistent  basis  sets,  are 
only  feasible  for  smaller  molecules  displaying  unusual  electronic  properties. 


19  This  calculation  was  run  with  the  standard  default  options  of  Gaussian98.  The  default  mode 
imposes  modest  convergence  criteria  sufficient  for  0.5  kj/ mol  accuracy  in  the  SCF  energy  and 
three  decimal  place  accuracy  in  quantities  derived  from  the  electron  densities  (e.g.  electrostatic 
potentials  and  derived  charges).  All  calculations  were  performed  under  the  IRIX  6.4  operating 
system  using  one  185MHz  CPU  optimised  for  floating-point  arithmetic. 

20  This  is  the  acronym  for  Dunning's  correlation-consistent  polarised  valence  triple-zeta  basis  set 
[67]  augmented  by  diffuse  functions.  'Correlation-consistent'  refers  to  optimisation  of  the  basis 
sets  with  electron-correlated  atomic  calculations. 

21  Recorded  CPU  time  for  single  RHF/aug-cc-pVTZ  energy  calculation  on  the  SGI  Origin  2000: 4 
days,  17  hours,  41  min  and  16.7  sec. 
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Density  functional  methods  avoid  the  computational  expense  of  other  post-HF 
methods  by  approximating  the  effects  of  electron  correlation.  DFT  calculation  times 
formally  scale  as  the  third  power  of  the  number  of  electrons  in  the  system,  Ne3,  whereas 
computational  times  for  MP4  or  QCISD(T)  scale  as  Ne5.  In  practical  terms,  this  means 
that  DFT  treatments  of  electron  correlation  in  TNAZ  should  only  take  four  times  as 
long  as  a  F1F  calculation22.  Given  its  demonstrated  accuracy  for  predicting  molecular 
energies  and  activation  barriers,  the  enormous  amount  of  time  saved  by  using  DFT 
makes  it  an  attractive  method  for  studying  the  ground  state  properties  and  reactions  of 
larger  energetic  molecules  like  RDX  or  HMX. 


4.3  High-accuracy/low-cost  model  chemistries 

Several  procedures  have  been  developed  to  accurately  predict23  thermochemical 
quantities  such  as  atomisation  energies,  ionisation  potentials,  electron  affinities  and 
proton  affinities,  at  relatively  low  computational  cost.  These  procedures  fall  into  two 
general  categories:  compound  methods,  and  complete  basis  set  (CBS)  models 

4.3.1  Compound  methods 

Compound  methods  such  as  Gaussian-1  [68]  and  Gaussian-2  [69]  ('Gl'  and  'G2') 
attempt  to  approximate  a  single  high-level  calculation  by  a  combination  of  several  low- 
level  calculations.  For  example,  the  procedure  for  a  Gl  calculation  is  as  follows: 

1.  Optimise  the  molecular  structure  using  HF/ 6-31G(d)  and  determine  the  vibrational 
zero-point  energy  (ZPE)24. 

2.  Reoptimise  the  molecular  structure  using  full  (all-electron)  MP2/  6-31G(d).  Perform 
all  subsequent  calculations  using  this  geometry. 

3.  Compute  a  base-level  energy  (Ebase)  at  MP4/ 6-311G(d,p). 

4.  Compute  the  energy  correction  AE+  obtained  by  including  diffuse  functions 
(MP4/6-311+G(d,p)):  AE+  =  EbaSe+diff  -  Ebase- 

5.  Compute  the  energy  correction  AE2df  obtained  by  including  higher  polarisation 
functions  on  atoms  with  Z  >2  (MP4/  6-311G(2df,p):  AE2df  =  EbaSe+2df  -  Ebase. 

6.  Compute  the  energy  correction  for  residual  electron  correlation  effects  (to 
counteract  the  known  deficiencies  of  truncating  perturbation  theory  to  fourth 
order)  using  QCISD(T)/6-311G(d,p):  AEQ0  =  Eqcisd(t)  -  EbaSe- 


22  Actual  test  time  for  single  B3LYP/6-31G(d)  energy  calculation  on  the  SGI  Origin  2000:  34  min 
18  sec. 

23  In  some  cases,  predictions  lie  within  1%  of  experimental  values. 

24  The  ZPE  is  generally  corrected  by  a  scale  factor  which  accounts  for  known  deficiencies  in  the 
harmonic  approximation  of  intramolecular  potential  energy  surfaces.  For  HF/6-31G(d),  the 
scale  factor  is  0.8929. 
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7.  Estimate  a  higher-level  correction  of  the  remaining  correlation  energy  between 
spin-paired  electrons  using  the  formula  AE™-0  =  -0.00019na  +  -0.00595np,  where  na 
is  the  number  of  spin-up  electrons  and  np  is  the  number  of  spin-down  electrons  in 
the  molecule. 

The  resulting  G1  energy,  EG1  =  Ebase  +  AE+  +  AE2df  +  AEQCi  +  AF.hic  +  ZPE,  is  essentially 
an  approximation  to  the  energy  of  a  much  more  computationally  intensive  calculation, 
QCISD(T)/6-311+G(2df,p). 

The  G2  procedure  includes  a  couple  of  corrections  to  the  G1  energy,  and,  for  molecules 
containing  first-  and  second-row  atoms,  is  reported  to  predict  properties  to  within  an 
average  absolute  deviation  of  6.7  kj/mol  of  experimental  results.  A  modified  version  of 
G2  known  as  G2(MP2)  uses  second-order  instead  of  fourth-order  perturbation  theory, 
providing  nearly  the  same  accuracy  as  full  G2  at  substantially  lower  computational 
cost  [1]. 

G1  and  G2  methods  have  been  used  to  study  the  smaller  energetic  molecules, 
nitromethane  (H3C-NO2)  [70,  71]  and  nitramide  (H2N-NO2)  [72],  For  nitromethane, 
heats  of  formation  and  gas-phase  acidities  calculated  with  G1  and  G2  show  very  good 
agreement  with  experiment  (see  Table  4.2),  despite  the  tendency  for  MP2  to 
overestimate  N-O  bonds  lengths.25  As  expected,  the  G2  predictions  are  slightly  better 
than  those  of  Gl. 

Table  4.2  Comparison  ofGl  and  G2  predictions  of  nitromethane  chemistry:  heats  of  formation 
and  gas  phase  acidities. 


AfH°298  (kj/mol) 

AandH0  (kj/mol) 

Gl 

-87.6 

1500 

G2 

-64.6 

1486 

_ 

-74.8  ±  1.0 

1491 ± 12 

4.3.2  Complete  basis  set  methods 

The  largest  errors  in  MO  thermochemical  calculations  result  from  the  deficiencies  of 
small  (finite)  basis  sets.  With  perturbative  methods  (MP2,  MP4,  etc.),  it  is  generally 
observed  that  successive  contributions  to  the  calculated  energies  decrease  with  the 
order  of  the  expansion,  so  that  the  sum  converges  to  an  asymptote  representing  the 
best  attainable  energy  for  a  particular  basis.  The  convergence  properties  are  easily 
approximated  if  the  MP2  correlation  energy  for  a  many-electron  system  is  rewritten  as 


25  Chemical  bonds  predicted  by  MP2  are  typically  too  long  because  MP2  exaggerates  the  loss  of 
electron  density  within  the  bonding  region.  In  an  extreme  case,  MP2/4-31G  calculations  of 
H2N-NO2  predict  a  N-N  bond  length  of  2.106  A,  as  compared  to  the  experimental  value  of  1.86 
A  [73], 
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a  sum  of  electron  pair  energies  in  a  particular  series  known  as  a  'pair  natural  orbital' 
(PNO)  expansion.  Complete  basis  set  (CBS)  methods  use  the  known  asymptotic 
convergence  of  PNO  expansions  to  extrapolate  the  energy  of  a  complete  basis  set26 
from  calculations  using  finite  basis  sets:  as  the  size  of  the  basis  set  increases,  the  level  of 
applied  theory  decreases.  Hence,  CBS  models  typically  include  a  HF  calculation  with  a 
very  large  basis  set,  an  MP2  calculation  with  a  smaller  basis  set,  and  one  or  more 
higher-level  calculations  with  a  medium-  to  modest-sized  basis. 

Two  typical  CBS  model  chemistries  are  CBS-4  and  CBS-Q.  CBS-4  yields  results  similar 
to  other  high-accuracy  methods  with  substantially  less  computational  cost;  however,  it 
is  noteworthy  that  the  optimised  geometry  is  obtained  using  a  rudimentary  calculation, 
HF/3-21G(d).  CBS-Q  generally  achieves  better  accuracies  than  G2,  and  is  also  less 
expensive  computationally.  CBS-4  and  CBS-Q  are  only  capable  of  handling  systems 
composed  of  first-  and  second  row  atoms.  More  detail  on  these  methods  is  given  in  [1]. 

4.3.3  Caveats  for  energetic  materials 

Care  should  be  taken  when  using  compound  methods  to  model  energetic  molecules 
containing  nitro  groups.  The  electronic  structures  of  nitrogen  oxides  display  strong 
correlation  effects  because  the  highest-occupied  and  lowest-unoccupied  molecular 
orbitals  (HOMO  and  LUMO)  have  similar  energies,  allowing  them  to  'mix.'  Thus,  an 
accurate  quantum  chemical  model  of  nitrogen  oxides  requires  not  only  a  description  of 
the  occupied  electron  orbitals  up  to  the  HOMO,  but  also  a  description  of  the  LUMO  - 
in  technical  jargon,  this  model  is  termed  a  'multiconfigurational'  treatment.  Any  model 
of  nitro  compounds  that  relies  predominantly  on  single-configuration  methods  (e.g., 
HF,  MP2,  MP4)  is  subject  to  error.  For  instance,  the  use  of  MP2  predictions  with  small 
basis  sets  grossly  overestimate  O-N  and  N-N  bond  distances  in  nitrogen  oxides,  and 
fail  to  predict  the  correct  dissociation  of  nitromethane  [74].  Such  errors  can  corrupt  the 
results  of  highly  parameterised  procedures  like  Gl,  G2  and  CBS-Q,  which  rely  on 
MP2/6-31G(d)  geometry  optimisations  [1].  These  errors  may  be  avoided  by  the  use  of 
other  procedures  such  as,  for  instance,  CBS-RAD  (complete  basis  set  approximation  for 
radicals),  which  uses  a  high-level  configuration-interaction  method,  QCISD/ 6-31G(d), 
to  determine  geometries  and  zero-point  energies,  followed  by  coupled-cluster  methods 
for  high-accuracy  calculations  [75].  To  date,  there  has  been  no  reported  use  of  CBS- 
RAD  for  studies  of  energetic  materials. 


5.  Conclusions 

Electronic  structure  calculations  provide  useful  estimates  of  the  energetic  properties  of 
chemical  systems,  including  preferred  molecular  structures,  spectroscopic  features  and 
probable  reaction  paths.  This  review  of  some  current  electronic  structure  calculations 


26  A  complete  basis  set  contains  an  arbitrary  number  of  terms,  and  can  therefore  generate  the 
'true'  molecular  wavefunction. 
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has  concentrated  on  discussion  of  ab  initio  techniques,  as  these  are  the  most  accurate 
and,  in  principle,  'universal'  methods.  Rapid  advances  in  computer  technology  are 
making  computationally  expensive  ab  initio  methods  increasingly  more  practical  for 
use  with  realistic  chemical  systems.  In  particular,  'cheaper'  methods  such  as  density 
functional  calculations  and  'layered'  models  such  as  ONIOM  are  continually  being 
refined,  and  show  promise  of  providing  consistent  and  accurate  chemical  predictions 
for  most  complex  systems.  Energetic  materials,  which  often  contain  azide  or  nitro 
groups  that  require  sophisticated  treatments  of  electronic  structure,  fall  into  this 
category.  A  subsequent  publication  [4]  will  describe  the  application  of  electronic 
structure  methods  to  estimates  of  the  chemical  properties  of  energetic  materials. 
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